pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse, stpp, skimr, gganimate, ggplot2, plotly, flexdashboard, DT,gridExtra, rlist, grid, animation)Take_Home Exercise 1: Geospatial Analytics for Social Good: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar
1.0 Introduction
The study of armed conflicts in Myanmar has gained critical importance in understanding the geographical distribution and intensity of violence across different regions. Myanmar’s complex ethnic composition and ongoing civil strife make it a unique case for geospatial analysis. This project aims to apply spatial and spatio-temporal point pattern analysis methods to uncover the patterns of armed conflict between January 2021 and June 2024.
By leveraging conflict data from the Armed Conflict Location & Event Data (ACLED) and geospatial tools, we will focus on visualizing and interpreting conflict density through heat maps, Kernel Density Estimation (KDE), and advanced spatio-temporal analysis.

Our analysis will focus on four types of conflict events:
- Battles,
- Explosion/Remote violence,
- Strategic developments,
- Violence against civilians,
with particular attention paid to quarterly patterns in conflict occurrence.
2.0 Setting Up The Environment & Dataset
2.1 Installing the required Packages
Key Packages Used in the Project:
sf: Handles simple features in R, allowing for spatial data manipulation and analysis. It is crucial for reading and managing geospatial data like shapefiles (e.g., Myanmar’s administrative boundaries).raster: Used for raster-based spatial data manipulation, especially for working with raster maps, such as Kernel Density Estimation (KDE) results.spatstat: A powerful package for spatial point pattern analysis. It helps to analyze and visualize spatial point data, particularly for identifying clusters or patterns in armed conflict events.sparr: Builds onspatstatand focuses on performing spatial and spatio-temporal kernel smoothing, which will be crucial for KDE and heatmap creation.tmap: A thematic mapping package that will allow us to create maps, including KDE visualizations on an OpenStreetMap base.tidyverse: A collection of data manipulation packages likedplyr,ggplot2, andpurrr. It’s essential for data cleaning, manipulation, and visualization tasks.stpp: Used for spatio-temporal point pattern analysis, crucial for analyzing how conflict events evolve in both space and time.skimr: A quick and comprehensive tool to provide summaries and descriptive statistics for datasets, helping in the initial exploration.gganimate: Extendsggplot2to create animated visualizations. We can use this for animated time-series or evolving conflict maps.ggplot2: The core plotting package in R, essential for creating visualizations like time series plots and KDE heatmaps.plotly: Useful for creating interactive visualizations, allowing users to explore spatial data interactively (e.g., hover over points to see conflict details).pacman: is a package management tool in R designed to streamline the process of loading and installing packages.
2.2 Data-set involved in this topic
For this analysis, we use two key datasets:
2.2.1 ACLED Armed Conflict Data
Location & Event Data (ACLED)platform, which maintains an extensive record of conflict events globally. For this specific analysis, we will limit the dataset by filtering based on the following parameters to streamline data preparation and minimize the need for extensive data cleaning:
| Data Parameter | Filter Category |
|---|---|
| Date Range | From 01/01/2021 to 30/06/2024. |
| Event Type | 1. Battles 2. Violence Against Civilians 3. Explosions/Remote Violence 4. Strategic Developments |
| Country | Myanmar |
ACLED Configuration Image

Code to Import ACLED Dataset
ACLEDData <- read_csv("data/raw/aspatial/2021-01-01-2024-06-30-Myanmar.csv")Rows: 42608 Columns: 35
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (20): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (15): year, time_precision, inter1, inter2, interaction, iso, latitude, ...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
2.2.1.1 Understanding the data set fields.
Referencing this ACLED Official codebook, this is the dataset that we are working with, not to bore you with the details are mainly interested in the following fields,
Event ID: Unique identifier for each conflict event.
Event Date: Date of occurrence.
Event Type: Type of conflict event (e.g., Battles, Remote Violence).
Latitude/Longitude: Coordinates of the event.
Fatalities: Number of fatalities resulting from the event.
Actors: The groups or individuals involved in the conflict (e.g., state actors, ethnic armed groups).
Admin Levels: Administrative region, district, and township where the event took place.
If you’re interested in the data set fields to explore more, here’s the full fields!
ACLED Full Table Fields Summary
| Fields name | Fields Description | Values |
| event_id_cnty | A unique alphanumeric event identifier by number and country acronym. This identifier remains constant even when the event details are updated. | E.g. ETH9766 |
| event_date | The date on which the event took place. Recorded as Year-Month-Day. | E.g. 2023-02-16 |
| year | The year in which the event took place. | E.g. 2018 |
| time_precision | A numeric code between 1 and 3 indicating the level of precision of the date recorded for the event. The higher the number, the lower the precision. | 1, 2, or 3; with 1 being the most precise. |
| disorder_type | The disorder category an event belongs to. | Political violence, Demonstrations, or Strategic developments. |
| event_type | The type of event; further specifies the nature of the event. | E.g. BattlesFor the full list of ACLED event types, see the ACLED Event Types table. |
| sub_event_type | A subcategory of the event type. | E.g. Armed clashFor the full list of ACLED sub-event types, see the ACLED Event Types table. |
| actor1 | One of two main actors involved in the event (does not necessarily indicate the aggressor). | E.g. Rioters (Papua New Guinea) |
| assoc_actor_1 | Actor(s) involved in the event alongside ‘Actor 1’ or actor designations that further identify ‘Actor 1’. | E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank. |
| inter1 | A numeric code between 0 and 8 indicating the type of ‘Actor 1’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). | 1, 2, 3, 4, 5, 6, 7, or 8. |
| actor2 | One of two main actors involved in the event (does not necessarily indicate the target or victim). | E.g. Civilians (Kenya)Can be blank. |
| assoc_actor_2 | Actor(s) involved in the event alongside ‘Actor 2’ or actor designation further identifying ‘Actor 2’. | E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank. |
| inter2 | A numeric code between 0 to 8 indicating the type of ‘Actor 2’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). | 0, 1, 2, 3, 4, 5, 6, 7, or 8. |
| interaction | A two-digit numeric code (combination of ‘Inter 1’ and ‘Inter 2’) indicating the two actor types interacting in the event (for more, see the section Actor Names, Types, and ‘Inter’ Codes). | E.g.3, 58 |
| civilian_targeting | This column indicates whether the event involved civilian targeting. | Either ‘Civilians targeted’ or blank. |
| iso | A unique three-digit numeric code assigned to each country or territory according to ISO 3166. | E.g. 231 for Ethiopia |
| region | The region of the world where the event took place. | E.g. Eastern Africa |
| country | The country or territory in which the event took place. | E.g. Ethiopia |
| admin1 | The largest sub-national administrative region in which the event took place. | E.g. Oromia |
| admin2 | The second largest sub-national administrative region in which the event took place. | E.g. ArsiCan be blank. |
| admin3 | The third largest sub-national administrative region in which the event took place. | E.g. MertiCan be blank. |
| location | The name of the location at which the event took place. | E.g. Abomsa |
| latitude | The latitude of the location in four decimal degrees notation (EPSG:32647). | E.g. 8.5907 |
| longitude | The longitude of the location in four decimal degrees notation (EPSG:32647). | E.g. 39.8588 |
| geo_precision | A numeric code between 1 and 3 indicating the level of certainty of the location recorded for the event. The higher the number, the lower the precision. | 1, 2, or 3; with 1 being the most precise. |
| source | The sources used to record the event. Separated by a semicolon. | E.g. Ansar Allah; Yemen Data Project |
| source_ scale | An indication of the geographic closeness of the used sources to the event (for more, see the section Source Scale). | E.g. Local partner-National |
| notes | A short description of the event. | E.g. On 16 February 2023, OLF-Shane abducted an unidentified number of civilians after stopping a vehicle in an area near Abomsa (Merti, Arsi, Oromia). The abductees were traveling from Adama to Abomsa, Arsi. |
| fatalities | The number of reported fatalities arising from an event. When there are conflicting reports, the most conservative estimate is recorded. | E.g. 3No information on fatalities is recorded as 0 reported fatalities. |
| tags | Additional structured information about the event. Separated by a semicolon. | E.g. women targeted: politicians; sexual violence |
| timestamp | An automatically generated Unix timestamp that represents the exact date and time an event was uploaded to the ACLED API. | E.g. 1676909320 |
2.2.2 Myanmar Administrative Boundaries (Shapefiles):
Obtained through Geonode Mimu, this shapefile helps us to build the map and set the boundary zone of each district of myanmar. This dataset provides the geographical boundaries of Myanmar’s administrative divisions, from the national level down to the township level. It is essential for mapping conflict events to specific regions.
Code to Import Shapefile Dataset
M_State_Sf <- st_read(dsn="data/raw/geospatial/stateLevel", layer = "mmr_polbnda_adm1_250k_mimu_1") Reading layer `mmr_polbnda_adm1_250k_mimu_1' from data source
`C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial\stateLevel'
using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84
M_State_SfSimple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84
First 10 features:
OBJECTID ST ST_PCODE ST_RG ST_MMR PCode_V
1 1 Ayeyarwady MMR017 Region ဧရာဝတီတိုင်းဒေသကြီး 9.4
2 2 Bago MMR111 Region ပဲခူးတိုင်းဒေသကြီး 9.4
3 4 Chin MMR004 State ချင်းပြည်နယ် 9.4
4 5 Kachin MMR001 State ကချင်ပြည်နယ် 9.4
5 6 Kayah MMR002 State ကယားပြည်နယ် 9.4
6 7 Kayin MMR003 State ကရင်ပြည်နယ် 9.4
7 8 Magway MMR009 Region မကွေးတိုင်းဒေသကြီး 9.4
8 9 Mandalay MMR010 Region မန္တလေးတိုင်းဒေသကြီး 9.4
9 10 Mon MMR011 State မွန်ပြည်နယ် 9.4
10 11 Nay Pyi Taw MMR018 Union Territory နေပြည်တော် 9.4
geometry
1 MULTIPOLYGON (((95.20798 15...
2 MULTIPOLYGON (((96.17964 19...
3 MULTIPOLYGON (((93.36931 24...
4 MULTIPOLYGON (((97.59674 28...
5 MULTIPOLYGON (((97.1759 19....
6 MULTIPOLYGON (((97.81508 16...
7 MULTIPOLYGON (((94.11699 22...
8 MULTIPOLYGON (((96.14023 23...
9 MULTIPOLYGON (((97.73689 15...
10 MULTIPOLYGON (((96.32013 20...
M_District_Sf <- st_read(dsn="data/raw/geospatial", layer = "mmr_polbnda_adm2_250k_mimu") Reading layer `mmr_polbnda_adm2_250k_mimu' from data source
`C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84
M_District_SfSimple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84
First 10 features:
OBJECTID ST ST_PCODE DT DT_PCODE DT_MMR PCode_V
1 1 Ayeyarwady MMR017 Hinthada MMR017D002 ဟင်္သာတခရိုင် 9.4
2 2 Ayeyarwady MMR017 Labutta MMR017D004 လပွတ္တာခရိုင် 9.4
3 3 Ayeyarwady MMR017 Maubin MMR017D005 မအူပင်ခရိုင် 9.4
4 4 Ayeyarwady MMR017 Myaungmya MMR017D003 မြောင်းမြခရိုင် 9.4
5 5 Ayeyarwady MMR017 Pathein MMR017D001 ပုသိမ်ခရိုင် 9.4
6 6 Ayeyarwady MMR017 Pyapon MMR017D006 ဖျာပုံခရိုင် 9.4
7 7 Bago (East) MMR007 Bago MMR007D001 ပဲခူးခရိုင် 9.4
8 8 Bago (East) MMR007 Taungoo MMR007D002 တောင်ငူခရိုင် 9.4
9 9 Bago (West) MMR008 Pyay MMR008D001 ပြည်ခရိုင် 9.4
10 10 Bago (West) MMR008 Thayarwady MMR008D002 သာယာဝတီခရိုင် 9.4
geometry
1 MULTIPOLYGON (((95.12637 18...
2 MULTIPOLYGON (((95.04462 15...
3 MULTIPOLYGON (((95.38231 17...
4 MULTIPOLYGON (((94.6942 16....
5 MULTIPOLYGON (((94.27572 15...
6 MULTIPOLYGON (((95.20798 15...
7 MULTIPOLYGON (((95.90674 18...
8 MULTIPOLYGON (((96.17964 19...
9 MULTIPOLYGON (((95.70458 19...
10 MULTIPOLYGON (((95.85173 18...
2.2.2.1 Understanding the data set fields.
| Field Name | Description |
| OBJECTID | This is a unique identifier for each feature in the dataset, typically used to identify individual records or polygons in the shapefile. |
| ST | This represents the State or Region in Myanmar. For example, in your dataset, “Ayeyarwady” refers to a state/region. |
| ST_PCODE | This stands for State Postal Code. It is a standardized code that represents each state or region in Myanmar, such as “MMR017” for Ayeyarwady. |
| DT | This stands for District or Township within the respective state/region. For example, “Hinthada” is a district or township within Ayeyarwady. |
| DT_PCODE | This stands for District/Township Postal Code. It is a standardized postal code for each district or township, such as “MMR017D002” for the Hinthada district/township in Ayeyarwady. |
| DT_MMR | This field could be the District/Township name in Myanmar script, written in the local language. It may be an alternative representation of the “DT” field, showing the name of the district/township in Myanmar’s native language. |
| PCODE_V | This could be a Version of the Postal Code or a verification value used internally in the dataset. In this case, the value is “9.4”, possibly indicating a specific version of postal codes or an accuracy measure. |
| geometry | This column represents the spatial data for each district/township. It contains the geometrical shape (MULTIPOLYGON) defining the boundaries of the state or district/township, with coordinates provided in longitude and latitude. |
3.0 Data Pre-processing
To ensure accuracy and usability of the data, several preprocessing steps will be undertaken for the different datasets.
3.1 Myanmar Shapefile
3.1.1 Setting the CRS for the
Since Myanmar uses CRS of 32647 and when we download the map it’s in WGS84, we should change it to 32647 .
st_crs(M_State_Sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_State_Sf <- st_transform(M_State_Sf, crs = 32647)
# Verify that the CRS has been correctly set
print(st_crs(M_State_Sf))Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]
st_crs(M_District_Sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_District_Sf <- st_transform(M_District_Sf, crs = 32647)
# Verify that the CRS has been correctly set
print(st_crs(M_District_Sf))Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]
3.1.2 Renaming and removal of column names
colnames(M_State_Sf) <- c("OBJECTID", "state","state_code","type", "state_myr", "mimi_version", "geometry")
M_State_Sf_Cleansed <- M_State_Sf %>% select(state, type, state_myr ,geometry)colnames(M_District_Sf) <- c("OBJECTID", "state", "state_code", "district", "district_code", "district_mmr", "mimi_version", "geometry")
M_District_Sf_Cleansed <- M_District_Sf %>% select(state, district, district_mmr, geometry)3.1.3 Checking for validity of data
When working with spatial data, it’s crucial to ensure that all geometries are valid. Invalid geometries can cause errors in analysis and visualization.
Checking Validity with
st_is_valid():Identifying Invalid Geometries:
Fixing Invalid Geometries with
st_make_valid()
#checking if it's valid
M_State_Sf_Validity <- st_is_valid(M_State_Sf_Cleansed)
M_State_Sf_Invalid <- which(!M_State_Sf_Validity)
if (length(M_State_Sf_Invalid) > 0) {
print("MPZ Invalid!")
print(M_State_Sf_cleansed[M_State_Sf_Invalid, ])
} else {
print("it's valid!")
}[1] "it's valid!"
#checking if it's valid
M_District_Sf_Validity <- st_is_valid(M_District_Sf_Cleansed)
M_District_Sf_Invalid <- which(!M_District_Sf_Validity)
if (length(M_District_Sf_Invalid) > 0) {
print("MPZ Invalid!")
print(M_District_Sf_cleansed[M_District_Sf_Invalid, ])
} else {
print("it's valid!")
}[1] "it's valid!"
3.1.4 Visualizing the mynamar map
On the top right, you can toggle between the district level and also the state level to understand more about the boundaries of Myanmar.
tm_shape(M_State_Sf_Cleansed) + # Base map (Myanmar boundaries)
tm_basemap(server = "OpenStreetMap.HOT") + # Add OpenStreetMap as the basemap
tm_polygons("state", # Color the base map by state
palette = "Set3", # Use Set3 color palette for states
border.col = "gray", # Border color for the states
alpha = 0.5, # Semi-transparent polygons
title = "State", # Legend title for states
legend.show = TRUE # Show legend for state colors
) +
tm_layout(main.title = "States of Myanmar", # Main map title
legend.outside = TRUE) # Position the legend outside the map
There is more than 80 district here, so it only showcases 30 :)
tm_shape(M_District_Sf_Cleansed) + # Base map (Myanmar boundaries)
tm_basemap(server = "OpenStreetMap.HOT") + # Add OpenStreetMap as the basemap
tm_polygons("district", # Color the base map by state
palette = "Set3", # Use Set3 color palette for states
border.col = "gray", # Border color for the states
alpha = 0.5, # Semi-transparent polygons
title = "District", # Legend title for states
legend.show = TRUE # Show legend for state colors
) +
tm_layout(main.title = "District of Myanmar", # Main map title
legend.outside = TRUE) # Position the legend outside the mapWarning: Number of levels of the variable "district" is 80, which is larger
than max.categories (which is 30), so levels are combined. Set
tmap_options(max.categories = 80) in the layer function to show all levels.
Some legend labels were too wide. These labels have been resized to 0.51, 0.48, 0.53, 0.52, 0.49. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

3.2 ACLED Data
3.2.1 Changing the Column Names
Since Myanmar’s regional hierarchy follows State, District, and Township levels, we will rename the columns accordingly:
admin1→ Stateadmin2→ Districtadmin3→ Township
This is important because different countries use different administrative hierarchies. For example, in Singapore, the hierarchy is organized by Region and Subzones. Adjusting these names ensures that our dataset aligns with Myanmar’s specific regional structure for accurate analysis.
ACLEDData_Cleanse <- ACLEDData %>%
select(event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1,
actor2, inter2, interaction, admin1, admin2, admin3, location, latitude,
longitude, fatalities) %>%
rename(state = admin1, district = admin2, township = admin3) %>%
mutate(across(where(is.character), ~ replace_na(.x, "NA")), # Replace NA in character columns with "NA"
across(where(is.numeric), ~ replace_na(.x, 0))) # Replace NA in numeric columns with 03.2.2 Adding a “Quarter-Year” Column
To facilitate our temporal analysis, we need to add a “quarter-year” column based on the event_date field. This can be done by adjusting the date format to represent the quarter and year, ensuring that each event is categorized by the specific quarter it occurred in (e.g., Q1-2021, Q2-2022). This will allow for easier analysis of conflict trends over time.
# Convert event_date to Date format (if it's not already a date)
ACLEDData_Cleanse$event_date <- as.Date(ACLEDData_Cleanse$event_date, format="%d-%b-%y") # Adjust the format if needed
# Add a new column that shows the quarter and year
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
mutate(quarter_year = paste0("Q", quarter(event_date), "-", year(event_date)))
head(ACLEDData_Cleanse)# A tibble: 6 × 18
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political viol… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
# ℹ 10 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
# district <chr>, township <chr>, location <chr>, latitude <dbl>,
# longitude <dbl>, fatalities <dbl>, quarter_year <chr>
3.2.3 Joining ACLED’s Codebook Description
ACLED’s stores their data for the column “interaction” and “inter1” and “inter2” in codes, using their code book, let’s reorganise their data for simplier view, we can reference the code book here to know what each code represent. Map it out as a csv file read it in and change accordingly.
3.2.3.1 Left joining inter1 and inter’s description.
For more details about each inter code read here.
ACLEDActorInterCode <- read_csv("data/raw/aspatial/ActorTypesInterCode.csv")Rows: 8 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
left_join(ACLEDActorInterCode, by = c("inter1" = "code")) %>%
rename(inter1_description = Description)
# Left join again for inter2
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
left_join(ACLEDActorInterCode, by = c("inter2" = "code")) %>%
rename(inter2_description = Description)
head(ACLEDData_Cleanse)# A tibble: 6 × 20
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political viol… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
# ℹ 12 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
# district <chr>, township <chr>, location <chr>, latitude <dbl>,
# longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
# inter1_description <chr>, inter2_description <chr>
3.2.3.2 Left joining interaction description.
For more details about each interaction code read here.
ACLEDInteractionCode <- read_csv("data/raw/aspatial/AcledInteractionCodes.csv")Rows: 44 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
left_join(ACLEDInteractionCode, by = c("interaction" = "code")) %>%
rename(interaction_description = Description)
head(ACLEDData_Cleanse)# A tibble: 6 × 21
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political viol… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
# district <chr>, township <chr>, location <chr>, latitude <dbl>,
# longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
# inter1_description <chr>, inter2_description <chr>,
# interaction_description <chr>
3.2.3 Making it a SF Object and reverse geolocate state and district
Since ACLED provides longitude and latitude data, I prefer to reverse geolocate the points to match Myanmar’s official state and district boundaries. We are uncertain how ACLED assigns these regions, so to ensure consistency, we remove the original state and district columns from the ACLED data and replace them with the geolocated values.
Steps:
Convert ACLED Data to an SF Object: Using longitude and latitude coordinates, transform the ACLED dataset into a spatial object. Remember thatt we have to set CRS 32647 here as well.
Perform a Spatial Join: Match the points from ACLED with the corresponding state and district boundaries from the
m_sfshapefile, selecting only those columns.Remove Original Columns: After the spatial join, drop the original
state,district, andtownshipcolumns from the ACLED dataset.Rename the Joined Columns: Rename the newly joined
state.yanddistrict.ytostateanddistrict, effectively replacing the original columns with the reverse-geolocated values.
# Step 1: Convert ACLEDDataCleanse to an sf object using longitude and latitude
ACLEDData_Cleanse_Sf <- st_as_sf(ACLEDData_Cleanse, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
# Step 2: Transform ACLEDDataCleanse_sf to the same CRS as m_sf (EPSG: 32647)
ACLEDData_Cleanse_Sf <- st_transform(ACLEDData_Cleanse_Sf, crs = 32647)
# Step 3: Perform a spatial join, selecting only the state and district from m_sf
reverse_geolocated_state <- st_join(ACLEDData_Cleanse_Sf, M_State_Sf_Cleansed[, c("state")], join = st_intersects)
# Step 2: Spatial join to add 'district' from M_District_Sf_Cleansed
reverse_geolocated_district <- st_join(reverse_geolocated_state, M_District_Sf_Cleansed[, c("district")], join = st_intersects)
# Step 3: Remove original 'state', 'district', and 'township' columns from ACLEDData_Cleanse (if they exist)
# This step removes the original columns, and then renames the newly joined columns
ACLEDData_Cleanse_Sf <- reverse_geolocated_district %>%
select(-contains("state.x"), -contains("district.x"), -contains("township")) %>% # Remove old state, district, township columns
rename(state = state.y, district = district.y) # Rename newly joined columns
ACLEDData_Cleanse <- st_drop_geometry(ACLEDData_Cleanse_Sf)
# View the updated data
print(ACLEDData_Cleanse_Sf)Simple feature collection with 42608 features and 20 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -208804.4 ymin: 1103500 xmax: 640934.5 ymax: 3042960
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 42,608 × 21
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political vio… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political vio… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political vio… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic dev… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political vio… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic dev… Strategic… Milit… 1 NA
7 MMR64325 2024-06-30 2024 Political vio… Battles Milit… 1 PSLF/…
8 MMR64326 2024-06-30 2024 Political vio… Battles PSLF/… 2 Milit…
9 MMR64328 2024-06-30 2024 Political vio… Battles Milit… 1 PSLF/…
10 MMR64330 2024-06-30 2024 Political vio… Battles Milit… 1 PSLF/…
# ℹ 42,598 more rows
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, location <chr>,
# latitude <dbl>, longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
# inter1_description <chr>, inter2_description <chr>,
# interaction_description <chr>, geometry <POINT [m]>, state <chr>,
# district <chr>
3.2.5 Visualizing it by Event Type
# Set tmap mode to plot for static maps
tmap_mode("view")tmap mode set to interactive viewing
# Create the tmap object with the base map and event markers
tm_basemap(server = "OpenStreetMap") + # Add OpenStreetMap as the base map
# Add Myanmar state boundaries with transparency
tm_shape(M_State_Sf_Cleansed) +
tm_polygons("state", alpha = 0.3, border.col = "gray",
title = "State Boundaries", legend.show = TRUE) + # Add a legend for state boundaries
# Add event markers (bubbles) with size based on fatalities
tm_shape(ACLEDData_Cleanse_Sf) +
tm_bubbles(size = "fatalities", # Marker size based on fatalities
col = "event_type", # Color markers by event type
palette = "Set1", # Use Set1 color palette for event types
border.col = "black", # Border color for bubbles
border.alpha = 0.5, # Semi-transparent border
title.size = "Number of Fatalities", # Title for bubble size legend
title.col = "Event Types", # Title for event type legend
legend.size.show = TRUE, # Show legend for bubble size
legend.col.show = TRUE) + # Show legend for event types
# Layout settings for the map, including title and legend positioning
tm_layout(main.title = "Myanmar's State Conflicts by Fatalities", # Main map title
main.title.size = 1.5, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.2, # Size of the legend title
legend.text.size = 1, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05)) # Inner margins for the maplegend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Legend for symbol sizes not available in view mode.